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Abstract 

In this paper we develop a convergence analysis in an infinite dimensional setting 
of the Levenberg-Marquardt iteration for the solution of a hybrid conductivity 
imaging problem. The problem consists in determining the spatially varying 
conductivity a from a series of measurements of power densities for various voltage 
inductions. Although this problem has been very well studied in the literature, 
convergence and regularizing properties of iterative algorithms in an infinite 
dimensional setting are still rudimentary. We provide a partial result under the 
assumptions that the derivative of the operator, mapping conductivities to power 
densities, is injective and the data is noise-free. Moreover, we implemented the 
Levenberg-Marquardt algorithm and tested it on simulated data. 

Keywords. Inverse problems, nonlinear ill-posed problems, iterative regular- 
ization, elliptic equations, hybrid imaging 
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1 Introduction 



A common problem in hybrid imaging consists in the determination of the 
spatially varying conductivity a > in a domain f2 C R™ from m measurements 
of power densities £i{<j) — er|Vuj(cr)| 2 (i — 1, . . . , m) inside f2 (resulting from 
m different injected currents fi). That is, the potentials U{ satisfy the elliptic 
equation 

div(frViti) = in f2 , (11) 

u i = fi on 9il. 

This problem is relevant, for example, in Acousto -Electrical Tomography (AET) 
|21[ I5| HB] and Impedance- Acoustic Tomography (IAT) [5]. 

We investigate the solution of this nonlinear inverse problem in an infinite- 
dimensional setting. We apply the Levenberg-Marquardt iteration [S], a well- 
known iteration method which we recap in the next section. Using recent results 
about the linearized power density operator [1], we analyze local convergence 
conditions of the iteration method (for sufficiently smooth a and noise-free data 
£ {) and provide a partial result. 

There are a number of theoretical results available on the problem of estimating 
er from power densities (and some additional boundary information). In a paper 
by Capdeboscq et al. [7] it was shown that in R 2 , the conductivity a is uniquely 
determined by measurements Ficr) = (cr|Vui| 2 , cr|Vit2| 2 , crVitiV^) if 



det(Vtti, Vu 2 ) > c > in n. (1.2) 

(certain Dirichlet boundary conditions are known to enforce this condition, see 
P]). Note that aVu{Wu2 can easily be obtained from a third measurement of the 
power density using the polarization identity. [5] extended this result to R 3 and 



|16j to arbitrary dimension (where the above determinant condition ( 1.2 ) is much 
harder to fulfil) and additionally showed Lipschitz-stability of the reconstruction 
for sufficiently regular conductivities. Using the same interior measurements 
as above, Lipschitz-stability of the linearized problem F'(a) : L 2 (Vt) ~ > L 2 (Q) 
(where a € C°°(f2) and C R 2 ) modulo the kernel of the linearized power 
density operator was shown by Kuchmcnt and Stcinhauer [151 . What may be 
reconstructed in the setting of only one measurement <x|Vw| 2 is analyzed in [5]. 
Numerically, the problem has been treated in j3J HH [T7] . 

The paper is organized as follows: First, we introduce the Levenberg-Marquardt 
iteration and its convergence conditions. Then, we linearize the power density 
operator £ with Dirichlet boundary conditions in a suitable topology, analyse its 
stability and injectivity and discuss convergence of the iteration method. In the 
last two sections, we describe our numerical implementation of the Levenberg- 
Marquardt method and present numerical results. 
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2 Iterative solution scheme 



2.1 The Levenberg-Marquardt iteration 

Let us denote by £ s = (£i)i=i,..., m the (noisy) measurements according to m 
different initializations. We want find a such that £(<r) ~ £ S . To solve, we use 
an iterative regularization method in a Hilbert space setting. In [TTJ U21 [TU] 
several such methods for solving a general operator equation F{a) ~ y s , where 
F: X — > Y and X, Y are Hilbert spaces, are given. We use a Newton- type 
method, since these methods are usually converging faster than pure gradient 
descent methods. Newton's method itself is defined by 

(Tk+i = °k + F'(a k y\y s - F(a k )) . 

If the operator F'(a k ) is not left-invertible (this can be the case for £(a) if there 
are not enough measurements, as we shall see in the next sections) one may use 
Tikhonov regularization: 

<7fc+i = argmhi - F(<j k ) - F' (a k )(a - CT fe )||y + a k \\a - cr k \\x\ ■ 
The resulting iteration method 

Tfc+i = cT k + {F'{cj k yF'{a k )+a k l&)- 1 F'{a k y{y s -F{a k )), (2.1) 

where a k > are regularization parameters, is called the Levenberg-Marquardt 
method. 

Hanke |9 suggests choosing a k as the solutions of 

\\y s - F{a k ) - F'(a k )(a k+1 (a k ) - a k )\\ y = q \\y s - F(a k ))\\ y (2.2) 

for some < q < 1. He shows that with the above choice of parameters and the 
initial value oo sufficiently close to the desired solution, the Levenberg-Marquardt 
method converges monotonically to a solution a s of F(a s ) = y s , provided the 
Frechet derivative F' is uniformly bounded in a ball B p (a s ) containing the initial 
value o-Q and 



\\ F{a) _ F{d) _ F > {a){a -<j)\\ Y <C\\<T-<,\\ x \\F(a) F{c)\\ Y (2.3) 
for all a, a e B p (a 5 ). 

2.2 Linearization 

We now prove Frechet differentiability of the power density operator £(a) = 



div(crVu) (where u solves (1.1 1 for some boundary condition / = /;). 
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Let C K n . We assume that a is positive and in the space H (fi) for I > 5 , which 
implies that H l (£l) is a Banach algebra with respect to pointwise multiplication 
PQ. For / G H l +^(dfl), (O) then has a unique solution u(cr) e H l+1 (fl) [E], 



so £: ff'(O) if'(fi). 

Let : fr — > H l ~ 2 , f 1-4 div(erV/). By formal differentiation one can see that 
the directional derivative u'(a)r is given by 

Lau'icr^T = — div(YV«) in f] 

(2.4) 

it (<t)t = on <9f2. 

Obviously, u'(a): H l (fl) —> H l+1 (£l) is linear with respect to r. It is also 
bounded for the given norms considering 



h'OHIffi+MO) < c 'll div ( T Vu)|| H i_i (n) < C||r|| Ki(a) . (2.5) 

Here we used the regularity estimates in [TI5] for the first and the Banach algebra 
property for the last inequality. Since the constant in the regularity estimates is 
a bounded function of a if a is bounded from below (e.g., in a sufficiently small 
ball around the solution), u'(a) is actually uniformly bounded there. Now, let 
i? u (cr, r) = u{<j + r) — u(cr) — m'(<t)t be the first order Taylor remainder of u. 
For small r we get 



L (7+T R u (a,T) = —L t u'{(j)t in fi. 



As in (2.5) we obtain 



\Ru{<y,T)\\ w+1(n) < C\\L T u'(a)T\\ H ^ 1 



< ciMlff'^) II u '( ct ) t IIh'+ i (^) - C W T 



i 2 



(2.6) 



since u'(a) is bounded. This shows that u'(a) as defined in (2.4) is actually the 
Frechet derivative of u: H l {tt) -> H l+1 (tt). 

For the power density operator £ : H l (£l) — > iJ'(Sl), a H> cr|VM(cr)| 2 , formal 
differentiation gives the directional derivative 



£'(a)T = |Vu| 2 t + 2ctVuVu'(<j)t, 

Its uniform boundedness follows from the uniform boundedness of u'(a) and the 
Banach algebra property of H l (Q). 

For the first order Taylor remainder Re (a, r), we get 
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R e (a, t) := £(a +r) - £{a) - £'{a)r 

= {<t + t)\V(u + u , {<t)t + R u (<j,t))\ 2 -a\Vu\ 2 
- \Vu\ 2 t - 2aVuVu'(a)r 

= (a + r)V(2tt + 2u'{g)t + R u (cr, r))VR u {a, t) 
+ (a + t)|Vu'(ct)t| 2 + 2tVmVm'(ct)t. 



Now, by (2.61 and (2.51 



<C(\\R T ) \\H l + 1 (n) + ll r Hif ! (n) + ll T ll// ! (o)) 
^ C ll T lli/i(n) 



(2.7) 



so \\Rs(a, t)\\ h i^ = odlrHjji/Qx), showing that £'(a) is the Frechet derivative 
off: H l (Q) -> if'(Q) at cr. 

2.3 Stability and injectivity 

In this section, we recap stability and injectivity of the linearized operator, 
following the treatment in jj] . 

We may write the linearization of S : H l (fl) — > H (ft; M. m ), a i-> (£i)i=i,..., m as 
the following system: 

V • 8aVu t + V • aVSui = 

(2.8) 

^cr|Viii| + 2crVwi • \7Sui = 5£i 

The last line provides the Frechet derivative S£i — d£i ■ 8 a seeing Sui = dui ■ 8 a 
as a function of 8a given by solving the first line with Dirichlct conditions. 

Ellipticity. Let us assume as before that a and Vui are in the space H l (fl) 
for I > ~. Then Dirichlet conditions for 8ui render the above system elliptic 
provided that qi(x,£) = for all 1 < j < m implies £ = 0, where 

\vui\ 

This theory is independent of dimension and is based on the theory of elliptic 
redundant systems of equations with boundary conditions satisfying the Lopatin- 
skii conditions; see [U [H5] . Similar ellipticity conditions are found in [T5] . From 
such a theory, we obtain for s — I with I > ^ , and coefficients a £ H s (n) and 
Ui G H s+1 (il), the following stability estimate 

||<5<7 - 5a\\ H ,(n) + \\Su - 5u\\ H s+i(niB.™) ^ ^ 

< C\\8£ - 8£\\ H s( n . Rm) + C 2 \\Su - 8u\\ L 2 (n . Rm y 
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Here, 8£ is the collection of terms 6 £i and Su the collection of Siii- The presence 
of the term with Ci indicates the fact that the systems may not be invertible. 
Rather, it is Fredholm, and hence invertible up to a finite dimensional kernel 
of (smooth) functions. This is reminiscent of the behavior of A + q(x), which 
is invertible up to a finite dimensional kernel, although it is invertible (with 
Dirichlet conditions) for most values of q(x). 



Injectivity of modified systems. Whether we can choose C2 = above is 
not known in general. Injectivity results are known in the presence of more 
measurements than necessary for ellipticity; see, e.g., [14] in the linearized setting 
and |3 rjj)] in the nonlinear setting. 

In the setting of two measurements in dimension n — 2 and three measurements 
in dimension n — 3, results of injectivity were obtained in pT for modified systems 



of equations. Let denote by v — {6a, {Sui}) and recast the system (2.8 1 as 



Av = S. 

Then assuming that full Cauchy data are known for v on dfl, we may replace 
the above system by 

A*Av = A*S in ft, v = d v v = on dfl. (2.10) 

This system was proved in [4J to admit a unique solution on sufficiently small 
domains; in other words, (|2.9|) holds with C2 = 0. 



A different modification of (2.8 1 also leading to a fourth-order system similar to 



(2.101 is also proved to be injective; see |U Equation (35)]. 



These systems are difficult to implement numerically. We therefore wish to solve 



the system as written in (2.8 1. Knowing that very similar systems are indeed 



shown to be injective in the aforementioned work, it seems reasonable to assume 



that C2 = in (2.9 1 for m large enough, i.e., that the linearized operator A is 
left-invertible. With this assumption, we get that d£ is an operator from H s (iY) 
to iF(ft; W m ) with a bounded left-inverse when s > |. 
We then obtain an optimal stability estimate of the form 

C- l \\Str - Sa\\ H . (a) < \\6£-5£\\ H s {u . nm) <C s \\5a-5a\\ H , m . (2.11) 
2.4 Convergence analysis 

In t his section, we analyse whether the convergence conditions given in section 



2.1 



can be shown to hold for the power density operator £ : H l (Q) — > iJ'(fi;K 2 ) 
with a suitable number of measurements m. Note that we require noise-free data 
here (the noise cannot be assumed fulfil differentiability constraints). 
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To get local convergence of the Levenberg-Marquardt iteration, it suffices to 
show that (2.3) holds (we already saw that the operator is uniformly bounded 
close to a solution). Given Lipschitz-type stability, i.e., 



H l (fl) 



<C'\\£(a 1 )-£(a 2 )\\ H , 



(2-12) 



the estimate (2.3 1 can easily be obtained from (2.7 1 



Using the stability estimates (2.11 1 for the linearized operator, we can get such 
an estimate locally, since for <7i,<72 € B p (cr s ) (a ball with sufficiently small 
diameter in the iJ'(f2)-norm) we have 

C ||<7 2 - < C7 1 \\<T2 - o"i||jj!(n) - c \\ a 2 - °"illj?'(n) 

< \\£'((ti)(<ti - o"2)|| H i (n . Km) -C\\a 2 - o-i||jj!(q) 

< 11^(0-2) - £ 0l)lliJl(f1;R*«) 

by using norm equivalence for the second inequality and 



\£'(ai)(<T 2 -<ri)\\ H i 



(fi;R m ) 



£(a 2 ) -^(o-i)|| Hi(a . H m) < C||o2-oi|| H1 



(n) 



which is obtained by applying the reverse triangle inequality to (2.7 1, for the 
third estimate. 



From (2.7 1 we finally get 



\Re((r,T)\\ H > 



(Q:i 



.) < C \\r\\ Hl(n) \\£{a + t) - S{o)\\ Hl{n . R „ 



(2.13) 



for sufficiently small r. Applying the general results in [S] now gives us local 
convergence (in the noise-free case) of the Levenberg-Marquardt iteration applied 
to S: H l (Q) -4 H l (n;R m ). 



Note that for the stability estimate (2.111 we assumed injectivity of (|2.8[) (and 



thus £' (a)). For one measurement, this may not hold, for example, it is easily 



seen that if Q, is the unit circle in R , the boundary voltage f(x,y) 



by 



(where ( %) is a unit vector) and a = 1, the function / (with corresponding 
linearized potential (u' (a)f)(x,y) = —\{x 2 + y 2 ) + \) is in Ker£'(er); see also 
0. 

Generically, one may however expect injectivity to hold if enough data are 
available. 



3 Numerical solution 

The following sections contain some numerical work to demonstrate the feasibility 
of our ideas. For simplicity, we work in two dimensions, so we take Vt C K 2 . 
Furthermore, we now allow the presence of noise. As noise cannot be assumed 
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to be differentiable, we have to consider (in contrast to the previous section) the 
power density operator to map into L 2 (Vi). Since I — 2 is the smallest integer 
for which H l (£l) is a Banach algebra, we take £ : H 2 (£l) — > L 2 (il). 

3.1 Calculation of the adjoint 

For the iterative algorithm, we require an expression for the adjoint of £'(&), 
which we now derive. 

First, we consider the case m — 1 of a single measurement. Larger values of m 
are treated with an additional summation. Since £\<j) : H 2 (Vt) — > L 2 {Vt) is a 
bounded operator, we know that £'{&)* : L 2 (Q) — > H 2 (Vt) is bounded as well 
(see [2D]). 

First note that the adjoint operator of £'(<r) : H 2 (fl) — > L 2 (n) can be written as 

£'{a)*: L 2 (n) -> H 2 (n) 
£'(a)*z = i*£'(a)*z, 

where £ {a)* is the L 2 (f2)-adjoint of £{u) and i* the adjoint of the embedding 
i: H 2 (n) -> L 2 (fl). 

We start by finding an expression for £ (a)* on H 2 {VL) C L 2 (f2). Let V — 
V(a,u(<r)) with V: H 2 (fl) — > Hq(Q) be the linear operator defined by 

L a Vz = — div (zaVu(a)) in f2 

= in an. 

Note that za is also in H 2 (£l), so V is mapped into L 2 (il). 
For z £ H 2 we now get that for the second summand of £'(a) 

(aVuVu'(a)T,z) L2{n) = - {aVu' (<j)t,Wz) L2{n . R2) 
= (T,X7uX7Vz) L2(n) 

by partial integration and the definitions of V and u'(a). Since the first summand 
is self-adjoint, we conclude 

£'(ay:H 2 (n)^H 2 (n) 

£'{a)*z = i*{\Vu{a)\ 2 z + 2Vu(a)VVz) . 

As £'{&)* is bounded on L 2 (fi), Q 

can be continuously extended to w € 
L 2 (Q), e.g., by taking £'(a)*w := lim e ^ £ '(cr)*$ e u> where $ e w = 4> e * w is the 
mollification operator corresponding to some mollifier 4>. 

For £: H 2 (£l) -4 L 2 (fi;M m ), a i-> (£i);=i m we similarly get (with Vj = 
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£'{a)*: i? 2 (0;K 2 ) -> iJ 2 (ft) 

m (3 2) 

S'(a)*z = i*Y^ l v %0)| 2 ^ + 2VujO) • VV>j . 1 ' ' 

3 = 1 



£' can also be continuously extended to L 2 (il; R 2 ) by mollifying and taking the 
limit. 

To calculate i* : L 2 (il) — > H 2 (Q.) we make the additional assumption that 
a £ H%(Cl) = {x e if 2 (f2)|§^ = 0}, a closed subspace of H 2 (Q) since o 
trace) : iT 2 (ft) -4 L 2 (<9ft) is bounded. 

With the inner product 

for some j3 > 0, i*y is the solution of the Neumann problem 

(ld+/3 2 AA)i*y = y 
dAi*y 

-^T = ° n ^ (3-3) 
ov 

since for all x £ H^(il), y £ L 2 (il) we get that 

{x,i*y) H 2 N(n) = (x,i*y) L2m + /3 2 (Ax,Ai*y) L2{n) 

= <», (W +/3 2 AA) l * y ) L2(0) + P 2 ^-Aey - ^JX 
= ( x ,y)L 2 (n) = i^'V) L 2 {n) ■ 



The fourth order PDE (3.3) is equivalent to the system of second order equations 



i*y + /3 2 Az = y 
Ai*y — z = 



£ = on 3tt ( 3 - 4 ) 
^ = 0onSa 

For a cruder approximation (which should give a faster solution), we could 
use H 1 or L 2 adjoints (where we have i = i* = Id). For the inner product 
( x >v)HHn) = 2/) L 2 (n) +/3(Va;,Vy) i2(0) , the adjoint of I: iT^Q) -> L 2 (f2; R 2 ) 
(by a proof similar to the above) maps y € L 2 (f2) to the solution of 
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(ld-f3 2 A)l*y = y 



di*y 
dv 



on an. 



(3.5) 



3.2 Implementation 



For F — £ with a given (single) measurement £ , the equation (2.1 1 for the fc-th 
Levenberg-Marquardt step r k reads 

(lim£*M e + a k lA) Tk = £'(a k )*(£ 6 - £ (a k )) (3.6) 

where, with u = u(a k ), a k = a, r = r k and a = o^., 

M £ t = |Vu| 2 $ £ (|Vu| 2 r + 2ctVwVm'(ct)t) 

+ 2VuVy$ e (r|Vu| 2 ) + 4VuVy$ e (crVMVu'(CT)T). 

Now, let e be fixed. By introducing auxiliary variables z 1 ,^ 2 ,^ 3 and j/ 1 , j/ 2 , j/ 3 



(for expressions which are a solution of a PDE) and using (3.4), equation (3.6) 
can be written as 

z 3 + are = £'(<?)* (£ s - 8(a)) 



Az a - z 2 = 0, 



/3 2 Az 2 + z 3 - z 1 =0, 
| Vw| 2 $ e (| Vm| 2 t £ + 20-VwVj/ 1 ) + 2VuVy 2 + 4VuVy 3 



dz 3 /^ = on 9^ 
dz 2 /dv = on 5f2 







(3.7) 



div($ e (a-VMVy 1 )CTVu) + L CT y 3 = 0, 
div($ e (T e |Vu|VVu) + L a y 2 = 0, 
div(r e Vu) + Lo-y 1 = 0, 



V = 



o on an 
o on an 
o on an 



We then obtain r k in (3.6) from r e by letting e — > (for numerical purposes the 
continuous extension, i.e., the mollification and the limit can be ignored, as it 
leaves the discretization unchanged in the limit e — > 0). Since its discretization 
does not contain inverse matrices (and thus can be solved using sparse matrix 
operations only), the system (3.7) is more amenable to numerical treatment 
than discretizing (3.6 1 directly. When using L 2 or H 1 adjoints, similar (smaller) 
systems can be generated. 

For multiple measurements £ , using the exact same approach leads to the system 
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2 3 + ar E =£'(a)*(S S -£{a)) 



Az 3 - z 2 = 0, dz 3 /dv = on dfl 

/3 2 Az 2 + z 3 - z 1 = 0, <9z 2 /di/ = on dQ. 



(3.8) 



|V«i| 2 $ e (|VM 3 | 2 T 6 + 20-ViijVyj 1 ) + 2VujVy 2 + 4VujVj$ - z 1 = 

3=1 

div($ e (crVM J Vy J 1 )crVu ; ,) + L ff y| = 0, y* = on <9fi 

div(<E> e (r e |V^| 2 )crVM J ) + L a y 2 = 0, y 2 = on <9fi 

div(T e Vw 3 -) + Lay] = 0, yj = on <9f2 

so 3 equations (and auxiliary variables) have to be added per additional mea- 
surement. We again obtain Tk from r £ by letting e — > 0. 



4 Results 



In our numerical solution we directly implemented the systems (3.7) and (3.8 1 to 
obtain Levenberg-Marquardt steps. For discretization of the partial differential 
equations we used a self-written linear finite elements framework in MATLAB. 
The triangular mesh (with 42849 nodes and 85007 elements) was created using 
DistMesh. 18J All calculations were done on a workstation computer. 

To test the reconstruction algorithm, we generated simulated data £ s = {£\, £\, £3) 
on the unit circle with the boundary conditions fi(x, y) — x, f%(x, y) — y and 
fs(x,y) = -^(x — y) (and added Gaussian noise with standard deviation 1 to 
avoid an inverse crime). 
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Conductivity a 



Power density £\ 





Power density £\ 



Power density £% 



Fig. 1: Conductivity and simulated measurements with boundary conditions 
fi(x,y) = x, f2(x,y) = y and fz{x,y) = A^(x — y) and Gaussian noise. The 
color axis was manually set to 1-10, the range of a. 
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For reconstruction, we set ctq = 1, /3 = 10 and chose ctk = (with a minimum 
of 10~ 8 ) a priori. Choosing au according to the criterion (2.2 1 would be possible 
(e.g., by reducing ak until the condition is met), but numerically expensive. 
Note that for two or more measurements, regularization is technically is not 
necessary since the operator can be assumed to be invertible (the linearized and 
discretized equations were solvable with ak = in all instances we tested), but 
advantageous for the iteration scheme since it helps ensure that iterates remain 
in a trust region and thus feasible (i.e., positive). Additionally, we enforce a 
minimal conductivity of 10 -12 on the iterates. 

We used the standard MATLAB sparse solver mldivide to solve the linearized 
and discretized equations. 

The results of our numerical experiments seen in Figures [2] and [3] show that from 2 
measurements, even in the presence of significant noise, very good reconstructions 
are possible. The third measurement, which considerably increases the runtime, 
mostly serves to slightly reduce the noise. The _ff 1 -approximation of the H 2 
adjoint works well without loss of accuracy. The L 2 -approximation, on the other 
hand, does not converge properly. From the difference images, we can see that 
in the main remaining error in the H 1 or H 2 reconstructions with 2 or more 
measurements is due to smoothing (which can be alleviated by running more 
iterations or lowering f3). 

While the given algorithms could in theory be directly translated to M 3 (a 
higher degree of regularity would be necessary to get a Banach algebra and/or 
map into L 2 (fl)), the drastically increased number of elements would make 
the computational effort for high-resolution reconstructions unreasonable. To 
reduce the effort, one could either switch to Landweber-like methods or use a 
combination of nested grids and iterative linear solvers as in [7j. 
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L 2 adjoints, data £ 2 



L 2 adjoints, data (£f,£|) L 2 adjoints, data (£\,£ s 2 ,£ s z ) 




H 1 adjoints, data £ 2 H 1 adjoints, data (£f,£g) H 1 adjoints, data (£f,£ 2 ,£^) 




H 2 adjoints, data £ 2 H 2 adjoints, data H 2 adjoints, data (£\,£ 5 2 ,£^) 

Fig. 2: Reconstructed conductivities after 15 iterations of the Levenberg- 
Marquardt algorithm using different adjoints and one, two or three measurements. 
The color axis was manually set to 1-10, the range of the original conductivity 
to ease comparison. 
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L 2 adjoints, data L 2 adjoints, data (£j,£|) L 2 adjoints, data (£^,£31 £3) 




H 1 adjoints, data S\ H 1 adjoints, data -H" 1 adjoints, data (£{,£\,£\) 




H 2 adjoints, data £^ H 2 adjoints, data (£^,£|) H 2 adjoints, data (£T|,£|,£|) 

Fig. 3: Difference images |<7 0S t — <r| corresponding to the reconstructions in Figure 
[2l The color axis was manually set to 0-3. 
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